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The continuum model of the twisted graphene bilayei^ is extended to include all types of com- 
mensurate structures. The essential ingredient of the model, the Fourier components of the spatially 
modulated hopping amplitudes, can be calculated analytically, for any type of commensurate struc- 
tures in the low twist angle limit. We show that the Fourier components that could give rise to a 
gap in the SE-even structures discussed by Mele- vanish linearly with angle, whereas the amplitudes 
that saturate to finite values, as — > 0, ensure that all low angle structures share essentially the 
same physics. We extend our previous calculations beyond the validity of perturbation theory, to 
discuss the disappearance of Dirac cone structure at angles below 6 < 1-. 
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I. INTRODUCTION 

Barely a year after the discovery of a new form of 
quantization of the Hall effect in graphene mono-layers 
layers^"—, the bilayer attracted considerable attention by 
displaying yet another type of Quantum Hall effect^. 
Experimental and theoretical studies quickly followed, 
on the electronic structure^. Landau level spectrurr^, 
transporter—, disorder and interaction s ^^i^^ . 

These early studies focused on the AB stacked 
bilayer—. Unlike the mono-layer, in which carriers near 
the Fermi level behave like massless fermions, the AB 
stacked bilayer has quadratic dispersion near the Fermi 
level (for undoped samples). It is gapless, as the mono- 
layer, but only in the absence of a perpendicular electric 
field. An important feature of this system is the existence 
of a variable energy gap induced by an external electric 
field perpendicular to the layer a^^'"'^^ . 

The first experimental indications of the existence of 
rotational disorder in ultra-thin graphite films came from 
films grown on the 4H - SiC(OOOi) (Carbon side) of SiC 
crystalsi^; however, it had been known for years that in 
graphite crystals the top layer is often found rotated with 
respect to the underlying ones, giving rise long wave- 
length modulations of the STM signals, displaying as 
Moire patterns22r— . Few-layer graphene films grown by 
chemical vapor deposition methods^ir— often show ro- 
tations of successive graphene layers. It has also been 
possible to produce twisted bilayers using mechanically 
exfoliated samples^l. 

The electronic structure of the twisted bilayer was first 
considered by the authors^ in the context of a contin- 
uum, Dirac- Weyl equation, description of the two lay- 
ers, coupled by a spatially modulated hopping. The 



model predicted the persistence of linear dispersion, with 
well defined Dirac cones, like in the mono-layer, but 
with an angle dependent suppression of the Fermi ve- 
locity; it was also predicted that there would be no gap 
in the presence of a perpendicular electric field. These 
results were subsequently confirmed experimentally by 
Raman22 and Landau level spectroscopy^, and by band 
structure calculations-^''^^, although the earliest calcula- 
tions appeared to question the suppression of the Fermi 
velocity^i^. The most striking confirmation of the elec- 
tronic structure proposed in^ came from the observation, 
with scanning tunneling spectroscopy, of two low energy 
Van-Hove peaks in the density of states, with a strongly 
angle dependent energy difference; these were identified 
with the occurrence of two saddle points in the band 
struct ure2^. 

The continuum description was originally developed 
for a specific family of commensurate structures, dense in 
the low angle limit, in which the relative displacement of 
corresponding Dirac points in each layer, AK = — K, 
(K^ is obtained from K by a rotation of the twist angle 
between the layers) is not a reciprocal lattice vector of the 
Moire super-lattice; as a consequence there is no direct 
hopping matrix element between these two Dirac points. 
MeleS considered the commensurability conditions more 
generally, and pointed out the existence of another fam- 
ily of structures in which AK is reciprocal lattice vector 
of the Moire super-lattice. This matrix element between 
the Dirac points of the two layers should then give rise to 
a significant gap, raising the possibility of quite different 
physics from the one discussed ini. 

Meanwhile, several authors22i2£i2^ addressed the 
physics at very low twist angles {9 < 1° ), finding sig- 
nificant deviations from some of the results presented in 
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our previous work. The continuum model is similar to a 
quasi-free electron calculation, where the kinetic energy 
scale is hvp^K = 2hvFKsm{9/2) and the periodic po- 
tential scale in given by the inter-layer hopping. The orig- 
inal calculation included a minimum set of plane waves, 
an approximation which in only valid if the kinetic energy 
scale dominates. 

In this work we review and extend the continuum 
model to address these issues. We are able to present 
a complete analytical calculation of all the Fourier com- 
ponents of the spatially modulated hopping for any fam- 
ily of commensurate structures in the low angle limit. 
The structures considered by Mele turn out to be quasi- 
periodic repetition of simpler structures of the type we 
originally considered. The Fourier components of the 
hopping amplitude that could lead to a gap, vanish as 
the angle decreases, due to an interference effect, whereas 
other amplitudes saturate, essentially ensuring that the 
low angle physics of all commensurate structures in the 
one we discussed previously. 

The complete characterization of the Fourier compo- 
nents of the interlayer hopping amplitude allows us to 
extend the treatment of the continuum model to very 
small angles. The Fermi velocity vanishes at an angle 
9 ^ 1- in very good agreement with the results obtained 
from band structure calculations22i22i^; an almost dis- 
persioneless band appears at this angle, corresponding to 
localized states around regions of AA stacking^. Using 
the continuum model, with only the dominant Fourier 
amplitude, Bistritzer and MacDonald'^'' showed that at 
even smaller angles the Fermi velocity becomes non-zero 
again, vanishing at a series of "magic angles", of which 
9 ^ 1- is the first in the series. We present a simple 
explanation of this observation based on the differences 
of the band structures of pure AB and pure AA stacked 
bilayers. 

In section |ll] we review the geometry of commensu- 
rate structures in the twisted bilayer in order to estab- 
lish notation and present a new derivation of the results 
obtained by Mele^ and Shallcross et. alr^. We formu- 
late the continuum model in section IIIII and present an 
analytical formulation of the calculation of the Fourier 
components of the spatially modulated inter-layer hop- 
ping, valid for small angles and any kind of structure. 
The main results of the model are presented in sectionlVl 
followed by a brief summary. 



II. GEOMETRY OF COMMENSURATE 
STRUCTURES 

The conditions for the commensurability of a Moire 
pattern of two rotated honeycomb lattices have already 
been considered by Mele^ and Shallcross et. al?'^ . We 
review this question, both to establish notation and to 
present an elementary approach to this question, more 
directly based on the symmetries of the hexagonal lat- 
tice. In this section we sketch the main argument, leaving 



details for Appendix 1X1 

The honeycomb (HC) lattice of graphene has an under- 
lying Bravais lattice with basis vectors which we choose 
as (lattice parameter a — 2.46 A) 



(la) 




(lb) 



This lattice is made up of two sub-lattices, A and B, 
where A atoms occupy Bravais lattice nodes, and the B 
are shifted by <5i = (ai + a2)/3: 



r^(m, n) = mai + na2 
rs(m, n) = ryi(m, n) + 5i 



(2a) 
(2b) 



In an AB stacked bilayer there are two such lattices, 
vertically displaced by c = 3.35 A, with the B atoms of 
layer 2 {B2) with the same horizontal positions as the A 
atoms in layer 1 (Ai), YB2{m,n) = YAi{'m,n). 

In a twisted bilayer, the layers are rotated relative to 
each other. We will assume we rotate layer 2 by an angle 
^, about a common A1B2 horizontal position, that we 
take to be the origin. A commensurate structure will 
occur if such a stacking A1B2 occurs elsewhere, say at 
Ti; the rotation might as well have been made about that 
second point, so Ti is a super-lattice translation, though 
not necessarily a primitive vector. For A1B2 stacking to 
occur, a B2 site must rotate to a Ai site, 

fcai + /a2 — > mai + na.2 k, I, m, n G Z, (3) 



which can only occur if 

fc^ + + fc/ = TO^ + + mn 



(4) 



since |fcai + /a2|^ ^ k"^ + P + kl. 

Shallcross et. al. in^^ present a detailed discussion 
of the solutions of this Diophantine equation. The same 
conclusions can be reached by exploring the point sym- 
metries of the hexagonal lattice, namely the existence of 
a six-fold rotation axis and of six reflection axis (the lines 
along the basis vectors ai,a2 and a2— ai, and three axis 
at angles of 7r/6 with these). These symmetries imply 
that a shell of Bravais lattice sites at a given distance 
from the origin, must be built of groups of two sets of 6 
sites, with position vectors. Pi and Q,;, « = 1, . . . , 6, such 
as displayed in Fig. ([Tb)) : the Pi(Qi) lie at directions 
making an angle of 7r/3, and the two sets are related to 
each other by reflection on the symmetry axis; these two 
sets may degenerate into one if it occurs on the symme- 
try axes. Naturally, a rotation of layer 2 by the angle 9 
that brings Pi — > Qi will leave six A1B2 sites at the Q 
sites, each defining a lattice translation of a commen- 
surate structure (from origin to Q^). The same can be 
said of the conjugate rotation 9' — Tr/3 ~ 9 that maps 
Qg — Pi, in which case the lattice translations are de- 
fined by the P^. Now, there may be, at a given shell. 
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first case, gcd(r, 3) = 1, the vertexes of the real-space 
Wigner-Seitz (WS) ceU of the super-lattice alternate be- 
tween B1A2 sites and hexagon centers; in the second case, 
each corner of the WS cell is an hexagon center of one 
layer and an atom of the other. In the reciprocal space, 
the shift in the Dirac point of the rotated layer, K^ — K, is 
a reciprocal lattice vector only in the second case. Mele^, 
who first called attention to these two types of commen- 
surate structures refers to them as sub-lattice exchange 
even (SE-even) when gcd(r, 3) = 3 and SE-odd when 
gcd(r,3) = 1. 



(a) 



(b) 



Figure 1: (a) Geometry of the Honeycomb lattice, (b) 
A shell of twelve Bravais lattice sites, their position 
related by the rotation and refiections symmetries of the 
hexagonal lattice. 



more than one of these groups of symmetry related sites. 
A shell of say 24 atoms will have two such groups P^, 
Qi and R^, S^. A rotation that, say, maps — > 
must map Si — >■ by symmetry, leaving us with 12 
A1B2 sites at the same distance from the origin: these 
lattice translations cannot be primitive translations, since 
the Bravais super-lattice is hexagonal by symmetry, and 
only has six nearest neighbors. Thus, in order to find all 
angles of commensuration, and the corresponding primi- 
tive vectors, we need only consider rotations that map 
{Pi} — ^ {Qi} or {Qi} {Pi}j where each of these 
sets of six points is obtained from the other by refiec- 
tion about the symmetry axes. 

These observations, and some elementary manipula- 
tions (see Appendix [5| are sufficient to establish the fol- 
lowing results for the possible commensurate structures. 

Angles: the following equation, with m and r co-prime 
positive integers, defines all possible angles of commen- 
surate structures with < 9 < n/3: 



cos 9{m, r) 



3™^ + 3mr + r'^/2 
3m2 + 3mr + r^ ' 



(5) 



Primitive vectors: the primitive vectors of the super- 
lattice for a commensurate structure of angle 9{m, r) are: 
i. Ifgcd(r,3) = 1, 



(6) 



tl' 




m 


m + r 




ai 


^2 




— (m + r) 


2m + r 




a2 



ii. If gcd(r, 3) = 3, 



"tl ■ 




.^2 





3 3 

■ m + 





ai 




a2 



(7) 



These two types of structures can be distinguished both 
in real and in reciprocal space-. Using the results of 
Appendix \^ it is straightforward to show that in the 



III. THE CONTINUUM MODEL 

The continuum description of the twisted bilayer was 
introduced by the authors^ in 2007. A single graphene 
layer admits an effective description in terms of the 
Dirac- Weyl equation for states close to one of the Dirac 
points^jiS. We use this description for the intra-layer 
Hamiltonians in the twisted bilayer, taking into account 
that layer 2 is rotated with respect to layer 1 by 9. We 
consider states near the Dirac point K = 47r(l,0)/3 in 
layer 1 and = (47r/3)(cos6', sin6') in layer 2. We de- 
note by ^'^(r), i — 1,2 the two component Dirac fields 
for each of the layers i = 1,2, and write the momentum 
as K -I- k in layer 1 and -I- k in layer 2. 

In momentum space the intra-layer Hamiltonians arei 

Hi -ft^*l k«i.r-k*i,k (8) 

k 

% -^^*2,k^VFT'-k*2,k; (9) 
fe 

the coordinate axes have been chosen to coincide with 
those of layer 1, r = (T^,Ty),T^ = e+»«^-/2Te-'^^-/2, and 
Tx and Ty are Pauli matrices. For the moment wc will ig- 
nore coupling between different Dirac valleys K, and 
K' = — K, K ^ = — ; we will return to this point 
later. 

To model the inter-layer coupling, Ji^ , we retain hop- 
ping from each site in layer 1 to the closest sites of layer 
2 in cither sub-lattice. We denote by "(r) the hori- 
zontal (in-plane) displacement from an atom of layer 1, 
sub-lattice a{a = Ai,Bi) and position r, to the clos- 
est atom in layer 2, sub-lattice f3' (/?' = A2,B'2)- The 
tight-binding inter-layer coupling is 

n^=J2t^ (^^'"(r,)) 4(rOc^, (r, + J^'"(r,)) + h.c. 

(10) 

where t± (^d°'^{r)j = t'^^{r), is the inter-layer, position 

dependent, hopping between orbitals with a relative 
displacement Co-|-<5, and Ca(j) is the destruction operator 
for the state in sub-lattice a at horizontal position r. 

Denoting by AK — — K the relative shift between 
corresponding Dirac wave vectors in the two layers, the 
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usual replacement Cq (r) 
to 



,1/2 



V'i,Q(r) exp(iK • r) leads 



a/3 

h.c. 



;K<'-5^"(r)g»AK.r,,,t 



V'I,a(r)V'2,/3 (r) 

(11) 



We used ?/'^(r + ^'^"(r)) « V'a(r) since the Dirac fields 
are slowly varying on the lattice scale. 

In Fourier space it is convenient to define ^i.k.a as the 
Fourier component of "02,0 (r) for momentum k± AK/2, 
the plus sign applying in layer 1 and the minus sign to 
layer 2. With this choice, the Dirac fields (/)i,k,a with 
the same k vector in both layers correspond to the same 
plane waves in the original lattice; the Dirac cones occur 
at k = -AK/2 in layer 1 and AK/2 in layer 2. 

For commensurate structures, the function 
t"''(r)exp 



is periodic and has nonzero 

Fourier components only at the vectors G of the 
reciprocal lattice: 



ifiG) 



1 

Vn 



(12) 



The integral is over the unit cell of the super-lattice, of 
area Vc- 

With these definitions the low energy effective Hamil- 
tonian, near K, is 



k,af3 



+ 



(G) 



2 J 
AK\ 



52, k,^ 



'2,k,, 



(13) 



Before proceeding it is perhaps worthwhile to remark 
that, including other interlayer hopping amplitudes, does 
not alter this description in a fundamental way. We 
would still arrive at a Hamiltonian similar to the one 
of Eq. Uni), but the hopping t'f (r) expliK'^ ■ 6°'^ {r) 

would be replaced by a more complicated expression. 

In this formulation, this problem is similar to that of 
a quasi-free electron band problem, because each layer 
has been reduced to a continuum, so that the only pe- 
riodicity remaining in the problem is that of the Moire 
super-lattice. The most important parameters are then 
the Fourier amplitudes, i^{G) defined by Eq. ([T^ . 

The implications of Mele's discussion of SE-even 
structures^ can now be clearly stated. In the SE-odd, r — 
1, structures we discussed in 2007, AK = (2Gi + G2)/3 
is not a reciprocal lattice vector of the Moire. There is 
no matrix element coupling between the Dirac cones K 
and K^ of the two layers. There is, in fact, a matrix el- 
ement coupling the different valleys, since K * — K is a 



G 





-Gi 


Gi G2 




i± 


i± 


i± 




i± 






ff^(G) 


i± 








i± 







Table I: The first and second line express exact results. 
In the next two lines these results have corrections of 
order a/L where L is the period of the super lattice; t± 
is real. 



reciprocal lattice vector; but this wave- vector has magni- 
tude 0(1/0), and for Moires with large periods, L ^ a, 
(r) is very slowly varying on the graphene lattice scale, 
and one would expect such matrix elements to be very 
small. But, as Mele pointed out, for an SE-even struc- 
ture, (gcd(r, 3) = 3), AK = r(Gi +G2)/3 is a reciprocal 
lattice vector of magnitude of order 0{1/L) and there 
seems to be no a priori reason to neglect it. It lifts the 
degeneracy between the two Dirac points and leads to a 
gap. A complete analysis of the Fourier amplitudes, to 
which we now turn, will allow us to resolve this issue. 



IV. CALCULATION OF FOURIER 
AMPLITUDES 

A. Structures with r = 1. 

We begin by considering the calculation of Fourier am- 
plitudes for r — 1 structures. Surprisingly, for small an- 
gles, the amplitudes for other structures can be reduced 
to these. 

In Reference^ we stated that in the low angle limit, 
and for an r = 1 structure, the dominant amplitudes are 
given by the results of Table HI We now give a complete 
justification of this statement, and show how one can 
calculate analytically all amplitudes for low angles. We 
begin by showing how certain symmetries imply relations 
between the horizontal shifts 6^^{r) for different sub- 
lattices. 

As stated in Section |lll three of the six vertexes of the 
WS cell are B1A2 sites: for instance, 



R 



2ti 



mai + Si ^ (m + l)a[ — S[. 



Since the origin is a A1B2 site, R is simultaneously a 
Ai — > Bi and a ^2 — > ^2 translation. Therefore, if there 
is an Ai site at r and B2 site at r -I- S^^{r), there will 
be to a Bi site at r -|- R and A2 site at r -|- R -|- 5 ' 
implying. 



r) + R = r + R + 5-''^(r + R), 



and ^^^(r) =^'^^(r-f R). 
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A somewhat more involved symmetry of this structure, 
namely, invariance under reflection about the origin, sub- 
lattice exchange {Ai o Bi, A2 o B2) and translation 
by R = msLi + Si, leads to a similar relation 5 (r) = 
—d^^{—r + R). These symmetries are exact and imply 
the following relations for the Fourier amplitudes: 



2BB (r-'\ -iG R (lAA 



tr(G) = e- 



(14a) 
(14b) 



With G = kGi + IG2, we get G • R = 27r(2/c - 0/3. 

The WS cell also has three vertexes which are hexagon 
centers (see Section|lI]); one such vertex is (ti + t2) /3 for 
r — 1 structures. This means that R = (ti + t2) /3 + 61 
is a Ai site and R' = (ti + 12) /3+<5i is an A2 site, and so 
5 AAi^ + Si) = 5'^ — 5i ~ 0{9). If this were exactly zero, 
R + ^1 would be a Ai Ai and B2 A2 translations, 
implying 



5AA{r) = 5BA{r - R - ^1) + 
This leads to 



i^^(G) « e-'^^-^^+'^ni^iG) 



(15) 



As before, G = fcGi + IG2 and G • (R + <5i) = 27r(fc + 
l)/i + 0{l/L). 

These three relations, Egs. p^ and (ITSI) . express all 
amplitudes in terms oii^^{G) and have been thoroughly 
confirmed by numerical evaluation of the Fourier ampli- 
tudes by calculating the integrals of Eq. as a lattice 
sum. For the specific values of G considered in Table |T] 
they lead to the phase factors relating amplitudes for dif- 
ferent sub-lattices. 

Let us now consider the expression for tf'^(G), and 
write it as a lattice sum. 



tr{G) 



-E 



In terms of G' : 
tf^(G) 



AK + G, we get 
1 



(16) 



ti_ [^^^(r.) 



^ giK''-(ri+<5B^(r,))g-»KTig-iG'Ti 



(17) 



In the WS cell of an r = 1 structure this simplifies be- 
cause exp iK^ • + S^^{ri)j exp [— iK • r] = 1 for all 

sites. It turns out, and this is a property that is exclu- 
sive to r = 1 structures, that the B2 site closest to Ai at 
= mai + na.2 is at tb' = msi'i + na.2 (same m and n), 
so that • tb' — K • r^. As a result. 



iG'- 



(18) 



Since tb^ ■= r^i + S^^{rA,) 



ta , where 



the rotation matrix, and for small angles, R{9) 



IS 

duix 



r, we get |5B^(r)| = 9r. We can therefore approximate 
the sum of Eq. [TS] as an integral 



tba 



-iG'-Ti 



Nr.n 



To simplify, we replace the hexagonal unit cell with a 
circle of the same area Nc<j = ^/3L'^/2, where L is the 
super-lattice parameter. The radius of the circle is i?ws = 

(V3/27r)^^^i and 



triG) 



V3L2 



dr 



' r cos 



47r 

71 ^ 



•^/3^ 



dxxtx^{x)jQ(G'Lx) (19) 



where Jq{x) is a Bessel function. 

To calculate this integral we need to parametrize the 
hopping between orbitals as a function of the horizon- 
tal shift 8. We express it in the Slater-Koster parameters, 
Vppa (d) and T^p^ (d) , where d is the distance between the 
two atomic centers, d = y^Cg + 6"^. For the d depen- 
dence of Vppcr{d) and YppT^^d) we used the parametriza- 
tion of ref.~; VppTj {ao / ^/3) , is the in-plane nearest neigh- 
bor hopping, t, and Vpp„{cQ) if the inter-layer hopping, 
t±, in an AB stacked bilayer. The contribution of I^ptt 
turns out to be negligible, and t±{d) is proportional to 
tj_: for 5 — ao/\/3, the carbon-carbon distance in a layer, 
<i((5)Ai«0.4. 

With this parametrization, we represent the amplitude 
as a function of G'L in Fig. [21 If t±{S) were constant, 
the integral would be proportional to Ji{G' L) / [C L) and 
decay as {G'L)~^^^. This is actually the way this ampli- 
tude decays, as could be seen by plotting G'^/^ times the 
integral. We have calculated numerically, as lattice sums, 
several amplitudes, using Eq. p^ : Fig.[2]shows that the 
analytical approximation to ^^^(G) gives an excellent 
account of the values found numerically. 

These results are worthy of the following comments: 

(i) The three reciprocal lattice vectors selected in Ta- 
ble |T1 G = 0, G = -Gi and G = -Gi - G2, ah have 
G'L = 47r/3. The corresponding values of i^^{G) are 
equal, f^'^(G) = OAt±; all other reciprocal lattice vectors 
have larger values of G", and the amplitudes are corre- 
spondingly smaller; these other amplitudes were ignored 
in Refs.i'^"^. 

(ii) For a general G = fcGi + ZGa, G' = (fc + 2/3)Gi + 
{I + 1/3)G2. Since G,; cx 1/L , G'L becomes independent 
of the angle or rotation. The amplitudes for a given {k, I) 
become independent of angle for small angles, tending to 
the values given by our analytical approximation. 

(iii) This complete characterization of the Fourier Am- 
plitudes, allows one, in principle, to include in the calcu- 
lation of the spectrum as many plane waves as necessary 
to achieve convergence. The characteristic energy from 
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1.0 



0.8 




3 3 3 3 3 3 

G'L 

Figure 2: The if ^(G)/ij_ as a function of G'L: the dots 
are numericany calculated values for a (m,r) = (10, 1) 
structure, with — 3.15°, and the red line is the 
integral of Eq. p^ . 



the in-plane motion is hvpAK ^ 0.190 9, with the energy 
in eV and the angle in degrees, and for small angles one 
requires more plane waves than those used in Ref.— . The 
physics of these small angle structures has been widely 
discussed recently in the literature and we will use these 
results to discuss it in the framework of the continuum 
model. But, before that, we consider the calculation of 
the Fourier amplitudes in other families of commensurate 
structures. 



B. Importance of r = 1 structures 

In this section we show that, in the small angle limit, 
the r = 1 structures are special, and determine the 
physics of all types of commensurate structures. 

In STM images^^. Moire patterns appear to satisfy the 
following relation between period and angle of rotation: 
L = a/ [2 sm{9/2)]. For a general (m, r) structure, 

/ 9(m, r)\ 1 r , ^ 

sin ^ ' ^ = - (20a) 
\ 2 J 2 Y 3m + 3mr + 

L(m,r) = a^/3m^ + 3mq + q^, (20b) 

where q = r/ gcd(r, 3), so the above relation is only satis- 
fied for r = 1. The plot 2Lsm{9/2)/a as a function 6, in 
Fig. |3al makes this clear. Remark that all these families 
of super-lattices, with different values of r, are dense as 
9^0. This means that a very small change in 9, with 
little effect in the structure in real space, can nevertheless 
change L by an arbitrary large factor. The implication 
is that, for very small angles, all commensurate struc- 
tures are almost periodic repetitions of structures with 
r = 1. That is seen very clearly by inspecting visually a 
few Moire patterns [see Fig. ((3b| ]. 

Let us show this explicitly for a SE-even structure, 
(m, r) with r = 3r'. At one of the corners of the Wigner- 
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(b) 



Figure 3: (a) 2{L/a) sin(6'/2) vs 9. The various lines 
correspond to different values of r; the lower line 
corresponds to the structures with r = 1. (b) A SE-even 
structure is almost periodic with the period of a 
structure with r — I; here is a (m, r) = (7, 3) is shown 
overlaid with the Wigner-Seitz Cells of (m, r) = (2, 1); 
the black hexagon is the true unit cell of the structure. 

Seitz cell, 

r — ^ — = "idi + -as = md 1 + -a^; (21) 

If TO mod 3 = 1, like in the (7,3) structure in Fig. I3bl 
this site has Bi atom of layer 1 and a hexagon center of 
layer 2. Therefore, at r — <5i there is an Ai site and at 
r — S'l, a B2 one. This implies that S^^{r) = S'l — Si — 
0{6). If this were zero, r would be a lattice translation of 
the Moire. The corresponding structure would be of SE- 
odd with m' = (to — l)/3 and r' = r/3. In real space, a 
SE-even structure (to, r), with to— 1 divisible by 3, is then 



very similar to a SE-odd with (m', r') = ((m— l)/3, r/3). 
In the following paragraphs we refer to these two lattices 
as £ (SE-even) and £ (SE-odd). 

Let us now relate the reciprocal lattice primitive vec- 
tors of £ and £ . Using the results of Appendix |X] one 
arrives at 



Gi 
G2 



2 1 

-1 1 



0(0) 



Gi 
G2 



(22) 



Ignore, for the moment, the 0{d) corrections. These 
equations tell us that the real space basis £, ti,t2, are 
linear combinations with integer coefficients of the basis 
of £. In the present case we have: 



(23) 



In the calculation of f^"(G) for the lattice with primitive 
vectors ti,t2, we can take into account that the 5^"(ri) 
are (approximately) periodic in ti and t2, and split the 
sum over in the unit cell £, into a sum over the unit cell 
of £, and a sum over the Uc unit cells of £ contained 
in the unit cell £: 



tl 




2 1 




tl 


_t2_ 




-1 1 




_t2_ 



3 

tl 



O -0.05 





1 




' 1 ' / 




0.4 


1 








0.2 








- 


















2 0.4 




- 










1 




,1 



tJG)/t^ (r = 1 structures, 2.5° < 9 < 7°) 

Figure 4: Comparison of Fourier amplitudes of pairs of 
structures: each point has an x coordinate ?^"(G) for 
((to — 1/3, r/3) and a y coordinate ^'^"(G) for a 
SE-even structure (m, r). G = kQ,\ + ZG2 and 
G = (2fc - /)Gi -I- (fc -I- OG2. According to Eq. ([26l) . 
these amplitudes should be equal. The line is y = x, not 

a fit. The inset has an expanded scale to include the 
dominant amplitudes (fc,Z) — {(0,0), (—1,0), (—1,-1)}. 
The angles are in the range 2.5° < 9 < 7.3°. 
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We achieved a factorization of t(["(G) 



(24) 



- 5(G) 



5(G) = 

Tic ^ 



-iG-T„ 



(25a) 
(25b) 



The second factor is t^"(G) for the lattice £. As for 
the structure factor, S'(G), note that, by definition, G • 
T = 2m7r, if T is a translation vector of £ (periodic 
boundary conditions), and, if G =G, a reciprocal vector 
of £, exp[iG • T„] — 1. Therefore we obtain, in this 
approximation. 



if(G) 



(26) 



where G is any reciprocal vector £. This is a very im- 
portant result: 

(i) : it expresses the Fourier amplitudes of SE-even in 

terms of those of structures with r = 1, which we 
calculated in section ITV Al 

(ii) : it states an approximate selection rule, that becomes 

more accurate as the angle of rotation decreases, al- 
lowing us to identify Fourier amplitudes that must 
tend to zero for small angles. 

We have checked this result by numerical calculation of 
<|["(G) for various lattices using Eq. ([TG]). In Fig.|4l each 
point on the plot has an x coordinate equal to t^"(G), 
G = A:Gi-|-;G2, for the {m\r') = ((m-l)/3, r/3) lattice, 
and a y coordinate i^["(G), G = (2/c - /)Gi + (fc + OG2 
of the SE-even (to, r) lattice; Eq. (|26p predicts that these 
amplitudes should be equal and the agreement is excel- 
lent. 

The second implication of Eq. (l26l) concerns the be- 
havior of Fourier Amplitudes for which G is not a re- 
ciprocal lattice vector of £, the r = 1 super-lattice. Of 
particular interest is G = AK = Gi + G2, because it 
determines the magnitude of the gap in a SE-even struc- 
ture. In Fig. [5] we show that f(["(AK) ^ as 6* ^ 0, 

as a result of the vanishing of the structure factor 5(G). 
In other words, there is a destructive interference in the 
sum of Eq. ([T6l) , because of the quasi-periodicity of the 
hopping amplitudes inside the unit cell of the larger pe- 
riod lattice. When G matches a reciprocal lattice vector 
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Figure 5: 



(G) /t. 



for G = kGi + IG2 for SE-even 

lattices, with different angles of rotation. For 
{k,l) = (—1, 1) or (2, 1), G is approximately equal to a 
reciprocal lattice vector of a r = 1 lattice and, 
i^{G)/t±\s almost constant; but for {k,l) = (1, 1), 
which corresponds to G = AK, the amplitude vanishes 
linearly with 6. 



G of the smaller period lattice S{G) « 1 (constructive 
interference) the amplitudes saturate to finite values as 
9^0. 

With this knowledge of the Fourier amplitudes for any 
structure, we finally address the calculation of the low 
energy bands of a small angle bilayer with a twist. 



V. THE CONTINUUM MODEL AT LOW 
ANGLES 

In the absence of the inter-layer coupling, 'H±, states 
with energy close to zero occur at k = — AK/2 in layer 1 
and k = +AK/2 in layer 2. The interlayer Hamiltonian 
'H± couples the states of momentum k in layer 1 to states 
k— G, in layer 2 with a matrix element (G) . The most 
important Fourier amplitudes, (of modulus i± = OAt±), 
in r = 1 structures, occur for G = 0, G = — Gi, and G = 
-Gi - G2 for which G'L = An/S, where G' = G + AK 
[see Fig. ([2])]. Neglecting other Fourier amplitudes^, the 
states of momentum k in layer 1 are coupled directly only 
to states of layer 2 of momentum k, k-|- Gi and k-l- Gi + 
G2 ; conversely the states of momentum k in layer 2 only 
couple to states k, k— Gi and k— Gi — G2. To investigate 
the spectrum at a momentum k close to zero energy, one 
can truncate the Hamiltonian to include only these six 
momentum values (three for each layer) giving a 12 x 12 
matrix to diagonalize (3 momentum values, 2 layers and 
2 sub-lattices)i. When k is close to the Dirac cone of 
one layer, the three momentum values that it couples to 
lie at the same distance AK from the Dirac point of the 
opposing layer; we have zero energy states coupling to 



two triplets of states at zLvpAK. The spectrum obtained 
from the diagonalization of the Hamiltonian matrix can 
be interpreted in a perturbative way when t^/vpAK <C 
1 . This analysis was presented in previous works and will 
not be repeated hereto. The main conclusions were: (i) 
the persistence of The Dirac cones, with linear dispersion; 
(ii) a renormalization of the Fermi velocity, relative to the 
single layer, which, in perturbation theory, was predicted 
as vp/vp = 1 — 9 (t±/{hvpAK)) {vp is the single layer 
value); (iii) the appearance of two low energy Van-Hove 
peaks due to the appearance of saddle points in the low 
energy bands, arising from the mixing of the two Dirac 
cones. 

In SE-even structures, however, there is a direct matrix 
element coupling the two Dirac cones; will the physics 
change relative to SE-odd structures due to the appear- 
ance of a gap? 

According to Eq. (1221) . the dominant Fourier ampli- 
tudes, in this case, occur for G = 0, G = — 2Gi — G2 
and G = — Gi — 2G2, since these correspond to G = 0,- 
Gi,-Gi - G2; on the other hand, AK = r(Gi + 
G2)/3 — Gi + G2. Therefore, these three dominant 
amplitudes couple the Dirac point of layer 1 to states 
of the layer 2 which are shifted from its Dirac point by 
-AK = -(Gi + G2) ,-AK + 2Gi -f G2 = Gi and 
— AK-t-Gi-|-2G2 = G2; since the angle between Gi and 
G2 is 27r/3, these are three vectors of the same modulus, 
|AK|, at 27r/3 angles, and we recognize exactly the same 
situation as discussed above for the r ~ 1 structures: 
the degeneracy points of each layer couple to two triplets 
at energies zLvpAK. It is true that, for this structure, 
there is a direct matrix element coupling the two degen- 
eracy points, corresponding to Gi + G2, which will lift 
the degeneracy and lead to a gap. However, as we saw in 
Fig. ([S]), this matrix element decreases with angle, and 
below 5° is under 5 meV. 

One can say that the differences between the various 
types of structures in momentum space are somewhat of 
a red herring. Two structures which are almost identical 
in real space must display similar physics. The momen- 
tum space description can look very different, but the 
magnitudes of the Fourier amplitudes must ensure simi- 
lar results. 

The perturbation theory in t±/ {hvpAK) clearly 
breaks down for very small angles, since, as we have seen 
the numerator becomes constant of order OAt± ~ 0.1 eV 
and the denominator is hvpAK ~ 0.190 x 6'eV (angle 
in degrees). This has led some authors22i22i^ to ques- 
tion the validity of the continuum description in the 
small angle limit. Band structure calculations, while 
confirming the prediction of depressed Fermi velocity, 
vp/vp = 1 — 9 (t±/{hvpAK)^ , find deviations from it 
below about 9 « 5°. 

One should not however confuse the perturbative re- 
sult with the continuum model. In fact, the continuum 
model should be better at smaller angles, since the scale 
of variation of the inter-layer hopping becomes larger. 
What one must do, however, is to include a larger set of 
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Figure 6: Velocity renormalization, by perturbation 
theory in t±/{hvF^K), and by the continuum model 
with numerical diagonalization of the secular equation, 

with a sufficient set of momentum values for 
convergence. The latter calculation deviates from the 
perturbative results as (hvpAK) becomes comparable 
to tj_, but is in quite good agreement with band 
structure calculationsSSi^iH. 



plane waves in order to achieve convergence of the low 
energy spectrum. 

In the following, we present some results for small an- 
gles, obtained by diagonalizing numerically the Hamilto- 
nian of Eq. p^ . truncated to a finite basis (largest ma- 
trix used of 168 x 168), and including all the required 
Fourier amplitudes, as given by the analytical expression 
of Eq. (fT9)) . This limit has already been addressed by 
Bistritzer and MacDonald^i in an approximation that 
includes only the dominant Fourier amplitudes. Some 
of our calculations, particularly those of the density of 
states, apparently require larger matrices for convergence 
than the ones that these authors claimed to have used. 

The results for the ratio of the Fermi velocity to the sin- 
gle layer value, as function of 6 are shown in Fig. ([5]). As 
expected, for small angles they deviate from the pertur- 
bative result, and compare very well with the values ob- 
tained from band structure calculations22i^2i^: the Fermi 
velocity becomes zero at about 9 ~ 1°. 

In Fig.[7a|we show a density plot of lowest positive en- 
ergy bands for 6 = 1.79° (vp/vp ~ 0.3); the Dirac cones, 
as well as the saddle point between them, are clearly visi- 
ble; at even smaller angles, = 1.20° , the corresponding 
plot shows an almost flat region in the arc joining the two 
Dirac cones through the saddle point (Fig. l7bl) : the range 
of energies with linear dispersion becomes very small. 

The density of states (DOS) is a very convenient tool to 
check for presence of Dirac cones in the band-structure. If 
the cones are present, the DOS shows a linear dependence 
near zero energy, as can be clearly seen in Fig. |5a] for 
= 1.79°; for — 1.2°, one can still define a (very small) 
Fermi velocity, but one should bear in mind that that the 
range of energies of linear dispersion is contracted to a 
few meV. 
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(a) (b) 

Figure 7: Contour plot of the first positive energy band; 
the hexagon is the First Brillouin Zone, (a) 9 — 1.79°, 

the cones are visible, but the saddle point is not located 
on the line joining the two Dirac cones; (b) — 1.20°, 
the arc joining the cones through the saddle point has 
become a very flat valley, and the cones are no longer 
well deflned. 
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Figure 8: Densities of states (DOS) for two angles; a) 
for 9 = 1.79° the cones are still well defined and the 
dispersion show the usual linear dependence near zero 
energy; b) for 9 = 1.2° there is a finite density of states 

near zero energy and one cannot define a Fermi 
velocity: the dispersion is no longer linear. The red line 
is the DOS for two uncoupled layers. 
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Figure 9: Densities of states (DOS) for two small 
angles; a) for 9 = 1.08° there is a very sharp central 

peak, with a barely resolved two peak structure, 
corresponding to a fiat band of states localized in AA 
stacking regions of the unit cell, b) for 9 = 0.87° the 
central band is broader, and still displays the two peak 
structure, although, at the meV resolution there is a 
finite DOS between the peaks, precluding the 
unambiguous definition of a Fermi velocity. 



At an even smaller angle, 9 = 1.08°, one observes a 
sharp peak in the DOS at low energy, corresponding to 
an almost dispersioneless band (Fig. [3a|), with a barely 
resolved two peak structure. Surprisingly, if the angle 
decreases further, the central band broadens (FIG: I9bp . 
This curious behavior was first found by Bistritzer and 
MacDonald^ and characterized as an oscillation of the 
Fermi velocity. In fact, at the meV resolution of the fig- 
ure, the DOS is finite between the peaks. It is not clear 
that a region of linear dispersion even exists, but, if it 
does, it is so narrow, that we prefer to concentrate on this 
curious variation of the width of the central peak. The 
diagonalization of the secular equation gives the eigen- 
states in the momentum basis, as well as the eigenvalues, 
so it is straightforward to calculate the local density of 
states. The density plots shown in Fig. [TU] show the local 
DOS in the superlattice unit cell, integrated over a nar- 
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Figure 10: Density plots of the Local Density of states, 
in the Wigner-Seitz Cell, integrated for |e| < 20 meV: 
(a) 9 = 1.78°;(b) 9 = 1.08°. 
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Figure 11: Band structures of pure A A and pure AB 
stacking bilayers. 



row energy range, close to zero; the dispersioneless band 
is composed of states localized in the AA stacking region, 
as was first found de Laissardiere and co-workers^. 

The reason for this localization, and for the curious 
fact that the degree of localization can oscillate with an- 
gle can be traced to the difference of band structures of 
A A and AB (or BA) bilayers. In a twisted bilayer of 
small angle, there are well defined regions of AA, AB 
and BA stacking, and it is legitimate to reason in terms 
of the corresponding band structures. The band struc- 
ture for an AA stacked bilayer, near each Dirac point, is 
composed of two cones shifted in energy by ±t±, corre- 
sponding to the bonding and anti-bonding combinations 
of Pz orbitals in each plane (Fig. [TT|) . As a result, the 
Fermi surfaces for electrons and holes, at zero energy, 
are circular with radius kp — t^/hvp- But on an AB 
or BA stacked bilayer at zero energy the Fermi surface 
is a point. From conservation of momentum parallel to 
an A A -H- AB interface, one can immediately conclude 
that there can be no transmission AA — AB{BA) for 
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any nonzero angle; in fact, a calculation shows that the 
transmission coefficient is also zero for zero angle of inci- 
dence and energy. Since the AA region is enclosed by an 
hexagon of AB and BA stacking, this raises the possibil- 
ity of localization of zero energy states in the AA region. 
However, this localization only occurs for zero energy; 
for finite energy some transmission is possible. Now, a 
confined AA region will have a discrete spectrum. The 
energy levels move toward the corresponding Dirac ener- 
gies (±ij_) if the size of the unit cell size increases. Is is 
clear then, that when a discrete level in an AA region oc- 
curs at zero energy, we can expect strong localization and 
a dispersioneless band. If we further decrease the angle, 
by increasing the unit cell, the discrete level moves away 
from zero energy, it starts tunneling into the neighboring 
AB and BA regions, and the low energy band broadens. 
This, we believe, is the rather simple explanation for the 
oscillation of the bandwidth of the central peak in the 
DOS. An explanation for this same oscillation, formu- 
lated in terms of non-Abelian effective gauge fields, has 
recently been proposed'^-. The extremely flat bands — 
which Bistritzer and MacDonalc^ associate with the ze- 
ros of the Fermi velocity — correspond to the passage 
of a state of a confined AA region through energies zLtj_ 
above (below) the corresponding Dirac energies {i.e. zero 
energy) . The double peak structures arises from the pres- 
ence of electron and hole states at zero energy in the AA 
regions. 



flat, band appears, separated by gaps from the rest of 
the spectrum (at positive and negative energies). This 
flat band is formed from states localized in AA stacking 
regions which, at zero energy, cannot tunnel into AB and 
BA regions. However, if the angle is further decreased 
the energy of these localized states changes, and they 
can start tunneling into the neighboring regions. This 
explains the oscillation with angle of the bandwidth of 
the central peak of the density of states. 
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VI. DISCUSSION AND CONCLUSIONS 

We have analyzed in detail the continuum descrip- 
tion of the twisted bilayer focusing on small angle struc- 
tures. We generalized our previous treatment to include 
all types of commensurate structures, and addressed in 
particular the possibility of a gapped electronic spectrum 
for SE-even structures raised by Mele^. We have shown, 
that for small angles, all commensurate structures are ei- 
ther of the type r = 1, in which the relation between the 
period and angle or rotation is that found in STM studies 
of Moire patterns, L — a/ [2sin(0/2)], or almost periodic 
repetitions of such structures. As a consequence, even 
though the momentum space description can be quite 
different, small angle commensurate structures share the 
same physics. 

We have achieved a complete analytical characteriza- 
tion of the Fourier components of the spatially modu- 
lated hopping amplitudes, which allows a detailed study 
of very small angle structures. This continuum descrip- 
tion accounts very well for the renormalization of the 
Fermi velocity relative to the single layer value. The den- 
sity of states is a revealing tool; if the angle is not two 
small, two well defined Van-Hove peaks appear at low 
energies, and, near zero energy, the DOS rises linearly, 
as expected for linear dispersion (Fig. [5a]); at 9 — 1.08- , 
the Van-Hove peaks are no longer resolved, as the range 
of linear dispersion shrinks to zero; a low energy, almost 



Given an arbitrary site of the hexagonal Bravais lat- 
tice. Pi — ka.1 + Isl2, the rotational/reflection symmetry 
implies that it is part of a set of twelve {Pi,Qi : i — 
1, . . . , 6}, 7r/3 being the angle between directions of con- 
secutive points in {P^} or in {Qi}, and each of these sets 
being the image of the other under reflection about the 
symmetry axes (see Fig. Ilbp . These two sets merge into 
one if and only if either k or I is zero or k = I. In the 
main text, we argued that we need only consider rotations 
that map one of these sets onto its image by reflection, in 
order to obtain all angles and primitive vectors of com- 
mensurate structures. Without loss of generality we can 
choose 

Pi = nsLi + msL2 (Ala) 
Qi = max + nsL2 (Alb) 
Qe = (to + n)ai — ma2. (Ale) 

with n > TO > 0; values of to or n zero, or m = n, cor- 
respond to 7r/3 rotations, that transform an AB stacked 
bilayer into an AA one (see also Fig. Ilb|) . So anticlock- 
wise commensurate rotations with angles < 6 < 7r/3, 
are of two types 

: Pi Qi; {n,m) {m,n); (A2a) 
6': Q'e^P'i; (p + g, -p) ^ (g,p); (A2b) 

In the flrst case Ti TOai +na.2 is a super-lattice trans- 
lation; in the second it is T' = qai +psi2 ■ We will soon see 
under what conditions these are primitive vectors. These 
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two rotations are conjugate, 6 + 9' = 7r/3, \i m = p and 
n = q. 

In the following, it will be useful to to define these 
rotations in terms of the pair of integers m, r with r = 
n — m, and p, s with s = q — p: 

(to + r, to) — > (to, to + r) (A3a) 
{2p + s,~p)^ {p + s,p). (A3b) 

One easily derives the following results for the angles, by 
taking the scalar product of final and initial vectors 



cos 6 



cos 9' 



Sto^ + Smr + 

Sto^ + 3mr + 
_ 3(TO + r/2)2 - (r/2)2 
^ 3(TO + r/2)2 + (r/2)2 

3pV2 + 3ps + 

3p2 + 3ps + 
_ 3(3p/2 + s)2 



(3p/2)2 



3(3p/2 + s)2 + (3p/2)2 



(A4a) 
(A4b) 
(A4c) 
(A4d) 



The second form of each expression makes it clear that 
the two families define the same set of angles: 9 — 9', 
\i m/r = s/Sp; all angles of commensurate structures 
are generated Eq. (jA4a[) with to and r positive integers: 
9'{p, s) = 9{m, r) if to = s and r = 3p. 

Given two positive integers, to, r, and the angle 9(m, r) 
defined by Eq. (|A4a|) . there is a unique set of integers p, q 
for which one of the following representations 



cos0(to, r) 
cos 9{m, r) 



ip^ + ?>pq + q^ /2 
3p2 + 2>pq + q^ ' 

3pV2 + 3pg + g2 
3p2 + "ipq + q2 



(A5a) 
(A5b) 



has the smallest denominator. If the smallest denomi- 
nator occurs for the first form, we conclude that ti :— 
p&i + (p + q)si2 is a lattice translation, with the smallest 
norm (the denominator is |ti| ) and, therefore, a prim- 
itive vector. The other can be obtained by a 7r/3 rota- 
tion of ti. On the other hand, if the second form has 
the smallest denominator, then, by the same reasoning, 
ti = (p + 'z)a-i + P3i.2 is a primitive vector of the super- 
lattice. 

From this point on, we assume that m, r are co-prime, 
because otherwise we can always reduce the denominator 
by factoring out the divisors of to and r. If 



3to2 



3TOr + rV2 _ ip^ + 2>pq + q^ 12 



(A6) 



3to2 -|- 3TOr -|- ip^ + ipq + q^ ' 

and 3^2 + 2tpq + q'^ < + 3TOr + r^, we must have, 

3to2 -I- 3TOr + r^/2 = X {3p^ + 3pq + q^/2) , (A7a) 
3to2 -I- 3TOr + = A (3^^ + 3pq + q^) , (A7b) 



where A is a positive integer. Subtracting these equa- 
tions, one gets = Xq^, so that X = s^, where s is a 



divisor of r. Solving the second equation for to/s, gives, 
recalling that m,r,p,q are positive integers. 



'"■ ^ _i_ ^ 



1 



+ 4p(p + q) ^p 



(A8) 



So s must a common divisor of to and r, and, since to, r 
are co-prime, s = 1, and the initial form already has 
the smallest denominator. An entirely similar argument 
can applied to reducing to the second form (Eq. (jA5b|) ). 
A form with smaller denominator is possible if r is a 
multiple of 3, and {p,q) — {m,r/3). 

In conclusion, we can state that if (to, r) are co-prime 
and 



COS! 



3to -|- 3TOr 



V2 



3to2 -f 3TOr + 



the super-lattice basis vectors are given by ti = SijSij, 
and the matrix S is defined in Eqs. ([6]) and ([7]). Shallcross 
et. al. define the angles and primitives vectors in terms of 
two co-prime integers p and g; their results coincide with 
these with the following correspondence: if r is odd, p = r 
and q = 2m -1- r; is r is even, p = r /2 and q — m ~\- r /2. 

From these results one can obtain other useful rela- 
tions. Since ti and t2 are lattice translations of both 
layers there have equally simple expressions in terms of 
the primitive vectors of the rotated layer sl^ and aj, . The 
transformation between these non-orthogonal basis is 



ai 




a2 





cos 9 + sin 9/^/3 -2sin6l/\/3 
2sin6'/\/3 cos 6* - sin 0/ %/3 



a'l 



(A9) 



The rotation matrix can be expressed in terms of to and 
r using Eq. (jA4ap . leading to 



tl 




m + r TO 




"a'l" 


t2 




—TO 2to -f- r 




a'2 



for 


gcd(r 


,3) = 


- 1 and 












tl 




TO 


+ 2r/3 


-r/3 




"a'l 






_t2_ 






r/3 


TO + r/3 






for 


gcd(r 


,3) = 


= 3. 











(AlO) 



(All) 



gl 


47r 


2 -1 




ai 




3|ai| 


-1 2 




a2 



The dual basis of {ai, a2} (reciprocal lattice primitive 
vectors) can be chosen as 



(A12) 



with a similar relation for {ti,t2} and its dual basis 
{Gi,G2}. Knowing that the Dirac points are given as 
K = (47r/3)(ai - a2), and K'^ = (47r/3)(a'i - a^, one 
can show, using Eqs. (TM)) to (|A12I) and Eqs. ((ST7)) . after 
some tedious but trivial algebra, the following relations. 



K" K 




if gcd(r, 3) = 1 
if gcd(r, 3) — 3. 

(A13) 

Note that, in the second case only, AK is a reciprocal 
lattice vector. 
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